#include "Matrice.h"

Matrice::Matrice ()  {

    for(int i=0; i<NB_LIGNES; i++){
        tab[i] = Vecteur();
    }
}

Matrice::Matrice (Vecteur liste[NB_LIGNES]){
    for(int i=0; i<NB_LIGNES; i++){
        tab[i] = liste[i];
    }
}


Matrice::Matrice(Vecteur v1, Vecteur v2, Vecteur v3, Point p)
{
	tab[0][0] = v1.x;
	tab[0][1] = v2.x;
	tab[0][2] = v3.x;
	tab[0][3] = p.x;

	
	tab[1][0] = v1.y;
	tab[1][1] = v2.y;
	tab[1][2] = v3.y;
	tab[1][3] = p.y;

	
	tab[2][0] = v1.z;
	tab[2][1] = v2.z;
	tab[2][2] = v3.z;
	tab[2][3] = p.z;

	
	tab[3][0] = v1.w;
	tab[3][1] = v2.w;
	tab[3][2] = v3.w;
	tab[3][3] = p.w;
}

int Matrice::getLignes(){return NB_LIGNES;}
int Matrice::getColonnes() { return NB_LIGNES;}

Point & Matrice::operator[] (unsigned i){
    return tab[i];
}

Point* Matrice::getTab(){ return tab;}

Matrice::~Matrice()
{
        
}

Matrice Matrice::multiplie(Matrice &m){

    Matrice prod;
	for(std::size_t lig=0;lig<NB_LIGNES;++lig){
		for(std::size_t col=0;col<NB_LIGNES;++col){
			//calcul du terme (lig,col)
			float tmp=0;
			for(std::size_t i=0;i<NB_LIGNES;++i){
				tmp += tab[lig][i]*m.tab[i][col];
			}
			prod[lig][col]=tmp;
		}
	}
	return prod;
}


Vecteur Matrice::multiplie (const Vecteur& v){

    return Vecteur(
				tab[0][0]*v.x + tab[0][1]*v.y+ tab[0][2]*v.z,// + tab[0][3]*v.w,
				tab[1][0]*v.x + tab[1][1]*v.y+ tab[1][2]*v.z,// + tab[1][3]*v.w,
				tab[2][0]*v.x + tab[2][1]*v.y+ tab[2][2]*v.z// + tab[2][3]*v.w
				);

}

Point Matrice::multiplie (const Point& v){

    return Point(
				tab[0][0]*v.x + tab[0][1]*v.y+ tab[0][2]*v.z + tab[0][3]*v.w,
				tab[1][0]*v.x + tab[1][1]*v.y+ tab[1][2]*v.z + tab[1][3]*v.w,
				tab[2][0]*v.x + tab[2][1]*v.y+ tab[2][2]*v.z + tab[2][3]*v.w
				);

}

Vecteur Matrice::translate(Vecteur original, Vecteur translation){

    Matrice m;
    m[0] = Vecteur(1, 0, 0);
	m[0][3] = translation[0];
    m[1] = Vecteur(0, 1, 0);
	m[0][3] = translation[1];
    m[2] = Vecteur(0, 0, 1);
	m[0][3] = translation[2];
    m[3] = Vecteur(0, 0, 0);
	m[0][3] = 1;

    return m.multiplie(original);

}
Matrice Matrice::matriceTranslate(Vecteur translation){

    Matrice m;
    m[0] = Vecteur(1, 0, 0);
	m[0][3] = translation.x;
    m[1] = Vecteur(0, 1, 0);
	m[1][3] = translation.y;
    m[2] = Vecteur(0, 0, 1);
	m[2][3] = translation.z;
    m[3] = Vecteur(0, 0, 0);
	m[3][3] = 1;

    return m;
}

	

Vecteur Matrice::rotation(Vecteur original, Point axe, float angle){

    Matrice m;
    m[0] = Point(
                    axe.x *axe.x*(1 - cos(angle))+cos(angle),
                    axe.x*axe.y*(1-cos(angle))-axe.z*sin(angle),
                    axe.x*axe.z*(1-cos(angle))+axe.y*sin(angle));

    m[1] = Point(
                    axe.y *axe.x*(1 - cos(angle))+axe.z*sin(angle),
                    axe.y*axe.y*(1-cos(angle))+cos(angle),
                    axe.y*axe.z*(1-cos(angle))+axe.x*sin(angle));

    m[2] = Point(
                    axe.x *axe.z*(1 - cos(angle))-axe.y*sin(angle),
                    axe.y*axe.z*(1-cos(angle))+axe.x*sin(angle),
					axe.z*axe.z*(1-cos(angle))+cos(angle));

    m[3] = Vecteur(0, 0, 0);
    return m.multiplie(original);

}

		

Matrice Matrice::matriceRotation(Vecteur const &axe, float angle){
	Vecteur vNormalise(axe);
	vNormalise.normer();
    Matrice m;
    m[0] = Vecteur(
                    vNormalise.x *vNormalise.x*(1 - cos(angle))+cos(angle),
                    vNormalise.x*vNormalise.y*(1-cos(angle))-vNormalise.z*sin(angle),
                    vNormalise.x*vNormalise.z*(1-cos(angle))+vNormalise.y*sin(angle));

    m[1] = Vecteur(
                    vNormalise.y *vNormalise.x*(1 - cos(angle))+vNormalise.z*sin(angle),
                    vNormalise.y*vNormalise.y*(1-cos(angle))+cos(angle),
                    vNormalise.y*vNormalise.z*(1-cos(angle))-vNormalise.x*sin(angle));

    m[2] = Vecteur(
                    vNormalise.x *vNormalise.z*(1 - cos(angle))-vNormalise.y*sin(angle),
                    vNormalise.y*vNormalise.z*(1-cos(angle))+vNormalise.x*sin(angle),
					vNormalise.z*vNormalise.z*(1-cos(angle))+cos(angle));

    m[3] = Point(0, 0, 0);

	return m;

}
	
Matrice Matrice::transpose(){
	Matrice m;
	int i,j;
	for(i=0; i<NB_LIGNES; i++){
		for(j=0; j<NB_LIGNES; j++){
			m.tab[i][j] = tab[j][i];
		}
	}
	return m;
}

Matrice Matrice::matriceRotation_x(float angle){

	Matrice m;
	m[0] = Vecteur(1,0,0);
	m[1] = Vecteur(0,cosf(angle),-sinf(angle));
	m[2] = Vecteur(0, sinf(angle), cosf(angle));
	m[3] = Point();

	return m;
}

Matrice Matrice::matriceRotation_y(float angle){

	Matrice m;
	m[0] = Vecteur(cosf(angle),0,sinf(angle));
	m[1] = Vecteur(0,1,0);
	m[2] = Vecteur(-sinf(angle), 0, cosf(angle));
	m[3] = Point();
	return m;
}


Matrice Matrice::matriceRotation_z(float angle){

	Matrice m;
	m[0] = Vecteur(cosf(angle), -sinf(angle),0);
	m[1] = Vecteur(sinf(angle),cosf(angle),0);
	m[2] = Vecteur(0, 0, 1);
	m[3] = Point();
	return m;
}

Matrice Matrice::operator*(Matrice& m){
	return multiplie(m);
}

Point Matrice::operator*(const Point& p){
	return multiplie(p);
}

Vecteur Matrice::operator*(const Vecteur& v){
	return multiplie(v);
}

Vecteur Matrice::changementRepere(Matrice repere, Point translation, Vecteur original){

    Matrice m;
    int i,j;
    for(i=0; i<3; i++){
        for(j=0; j< 3; j++)
            m.tab[i][j] = repere.tab[i][j];
        m.tab[i][j] = translation[i];
    }
    m.tab[i]= Vecteur(0,0,0);

    return m.multiplie(original);
}



void Matrice::afficher(std::ostream &flux) const {

    for(int i=0; i< NB_LIGNES; i++)
    {
           
            flux << tab[i].x<<" "<< tab[i].y<<" "<< tab[i].z<<" "<<tab[i].w<<" " ; //<<tab[i][1]<<" "<<tab[i][2]<<" "<<tab[i][3];
			flux<<std::endl;
    }
}

Matrice Matrice::matriceInverse()
{
	Matrice matriceInverse;


	float detA = tab[0][0]*tab[1][1]*tab[2][2]*tab[3][3] + tab[0][0]*tab[1][2]*tab[2][3]*tab[3][1] + tab[0][0]*tab[1][3]*tab[2][1]*tab[3][2]
				+tab[0][1]*tab[1][0]*tab[2][3]*tab[3][2] + tab[0][1]*tab[1][2]*tab[2][0]*tab[3][3] + tab[0][1]*tab[1][3]*tab[2][2]*tab[3][0]
				+tab[0][2]*tab[1][0]*tab[2][1]*tab[3][3] + tab[0][2]*tab[1][1]*tab[2][3]*tab[3][0] + tab[0][2]*tab[1][3]*tab[2][0]*tab[3][1]
				+tab[0][3]*tab[1][0]*tab[2][2]*tab[3][1] + tab[0][3]*tab[1][1]*tab[2][0]*tab[3][2] + tab[0][3]*tab[1][2]*tab[2][1]*tab[3][0]
				-tab[0][0]*tab[1][1]*tab[2][3]*tab[3][2] - tab[0][0]*tab[1][2]*tab[2][1]*tab[3][3] - tab[0][0]*tab[1][3]*tab[2][2]*tab[3][1]
				-tab[0][1]*tab[1][0]*tab[2][2]*tab[3][3] - tab[0][1]*tab[1][2]*tab[2][3]*tab[3][0] - tab[0][1]*tab[1][3]*tab[2][0]*tab[3][2]
				-tab[0][2]*tab[1][0]*tab[2][3]*tab[3][1] - tab[0][2]*tab[1][1]*tab[2][0]*tab[3][3] - tab[0][2]*tab[1][3]*tab[2][1]*tab[3][0]
				-tab[0][3]*tab[1][0]*tab[2][1]*tab[3][2] - tab[0][3]*tab[1][1]*tab[2][2]*tab[3][0] - tab[0][3]*tab[1][2]*tab[2][0]*tab[3][1];

	float idetA = 1 / detA;
	matriceInverse[0][0] = idetA * (tab[1][1]*tab[2][2]*tab[3][3] + tab[1][2]*tab[2][3]*tab[3][1] + tab[1][3]*tab[2][1]*tab[3][2] - tab[1][1]*tab[2][3]*tab[3][2] - tab[1][2]*tab[2][1]*tab[3][3] - tab[1][3]*tab[2][2]*tab[3][1]);
	matriceInverse[0][1] = idetA * (tab[0][1]*tab[2][3]*tab[3][2] + tab[0][2]*tab[2][1]*tab[3][3] + tab[0][3]*tab[2][2]*tab[3][1] - tab[0][1]*tab[2][2]*tab[3][3] - tab[0][2]*tab[2][3]*tab[3][1] - tab[0][3]*tab[2][1]*tab[3][2]);
	matriceInverse[0][2] = idetA * (tab[0][1]*tab[1][2]*tab[3][3] + tab[0][2]*tab[1][3]*tab[3][1] + tab[0][3]*tab[1][1]*tab[3][2] - tab[0][1]*tab[1][3]*tab[3][2] - tab[0][2]*tab[1][1]*tab[3][3] - tab[0][3]*tab[1][2]*tab[3][1]);
	matriceInverse[0][3] = idetA * (tab[0][1]*tab[1][3]*tab[2][2] + tab[0][2]*tab[1][1]*tab[2][3] + tab[0][3]*tab[1][2]*tab[2][1] - tab[0][1]*tab[1][2]*tab[2][3] - tab[0][2]*tab[1][3]*tab[2][1] - tab[0][3]*tab[1][1]*tab[2][2]);
 
	matriceInverse[1][0] = idetA * (tab[1][0]*tab[2][3]*tab[3][2] + tab[1][2]*tab[2][0]*tab[3][3] + tab[1][3]*tab[2][2]*tab[3][0] - tab[1][0]*tab[2][2]*tab[3][3] - tab[1][2]*tab[2][3]*tab[3][0] - tab[1][3]*tab[2][0]*tab[3][2]);
	matriceInverse[1][1] = idetA * (tab[0][0]*tab[2][2]*tab[3][3] + tab[0][2]*tab[2][3]*tab[3][0] + tab[0][3]*tab[2][0]*tab[3][2] - tab[0][0]*tab[2][3]*tab[3][2] - tab[0][2]*tab[2][0]*tab[3][3] - tab[0][3]*tab[2][2]*tab[3][0]);
	matriceInverse[1][2] = idetA * (tab[0][0]*tab[1][3]*tab[3][2] + tab[0][2]*tab[1][0]*tab[3][3] + tab[0][3]*tab[1][2]*tab[3][0] - tab[0][0]*tab[1][2]*tab[3][3] - tab[0][2]*tab[1][3]*tab[3][0] - tab[0][3]*tab[1][0]*tab[3][2]);
	matriceInverse[1][3] = idetA * (tab[0][0]*tab[1][2]*tab[2][3] + tab[0][2]*tab[1][3]*tab[2][0] + tab[0][3]*tab[1][0]*tab[2][2] - tab[0][0]*tab[1][3]*tab[2][2] - tab[0][2]*tab[1][0]*tab[2][3] - tab[0][3]*tab[1][2]*tab[2][0]);
 
	matriceInverse[2][0] = idetA * (tab[1][0]*tab[2][1]*tab[3][3] + tab[1][1]*tab[2][3]*tab[3][0] + tab[1][3]*tab[2][0]*tab[3][1] - tab[1][0]*tab[2][3]*tab[3][1] - tab[1][1]*tab[2][0]*tab[3][3] - tab[1][3]*tab[2][1]*tab[3][0]);
	matriceInverse[2][1] = idetA * (tab[0][0]*tab[2][3]*tab[3][1] + tab[0][1]*tab[2][0]*tab[3][3] + tab[0][3]*tab[2][1]*tab[3][0] - tab[0][0]*tab[2][1]*tab[3][3] - tab[0][1]*tab[2][3]*tab[3][0] - tab[0][3]*tab[2][0]*tab[3][1]);
	matriceInverse[2][2] = idetA * (tab[0][0]*tab[1][1]*tab[3][3] + tab[0][1]*tab[1][3]*tab[3][0] + tab[0][3]*tab[1][0]*tab[3][1] - tab[0][0]*tab[1][3]*tab[3][1] - tab[0][1]*tab[1][0]*tab[3][3] - tab[0][3]*tab[1][1]*tab[3][0]);
	matriceInverse[2][3] = idetA * (tab[0][0]*tab[1][3]*tab[2][1] + tab[0][1]*tab[1][0]*tab[2][3] + tab[0][3]*tab[1][1]*tab[2][0] - tab[0][0]*tab[1][1]*tab[2][3] - tab[0][1]*tab[1][3]*tab[2][0] - tab[0][3]*tab[1][0]*tab[2][1]);
 
	matriceInverse[3][0] = idetA * (tab[1][0]*tab[2][2]*tab[3][1] + tab[1][1]*tab[2][0]*tab[3][2] + tab[1][2]*tab[2][1]*tab[3][0] - tab[1][0]*tab[2][1]*tab[3][2] - tab[1][1]*tab[2][2]*tab[3][0] - tab[1][2]*tab[2][0]*tab[3][1]);
	matriceInverse[3][1] = idetA * (tab[0][0]*tab[2][1]*tab[3][2] + tab[0][1]*tab[2][2]*tab[3][0] + tab[0][2]*tab[2][0]*tab[3][1] - tab[0][0]*tab[2][2]*tab[3][1] - tab[0][1]*tab[2][0]*tab[3][2] - tab[0][2]*tab[2][1]*tab[3][0]);
	matriceInverse[3][2] = idetA * (tab[0][0]*tab[1][2]*tab[3][1] + tab[0][1]*tab[1][0]*tab[3][2] + tab[0][2]*tab[1][1]*tab[3][0] - tab[0][0]*tab[1][1]*tab[3][2] - tab[0][1]*tab[1][2]*tab[3][0] - tab[0][2]*tab[1][0]*tab[3][1]);
	matriceInverse[3][3] = idetA * (tab[0][0]*tab[1][1]*tab[2][2] + tab[0][1]*tab[1][2]*tab[2][0] + tab[0][2]*tab[1][0]*tab[2][1] - tab[0][0]*tab[1][2]*tab[2][1] - tab[0][1]*tab[1][0]*tab[2][2] - tab[0][2]*tab[1][1]*tab[2][0]);


	return matriceInverse;
}

Matrice Matrice::matriceHomotethie(Vecteur homotethie)
{
	Matrice matriceHomotethie;

	matriceHomotethie.tab[0][0] = homotethie.x;
	matriceHomotethie.tab[1][1] = homotethie.y;
	matriceHomotethie.tab[2][2] = homotethie.z;
	matriceHomotethie.tab[3][3] = 1;

	return matriceHomotethie;
}

Matrice Matrice::matriceRepere(Vecteur vecteurZ, Point &p)
{
	Vecteur vecteurX = Vecteur();
	float r1, r2, r3;

	do
	{
		r1 = (rand () * 2.0f - 1.0f) / RAND_MAX;
		r2 = (rand () * 2.0f - 1.0f) / RAND_MAX;
		r3 = (rand () * 2.0f - 1.0f) / RAND_MAX;
		vecteurX = Vecteur(r1, r2, r3);
	}while((vecteurX.pScalaire(vecteurZ)) == 0);
 
	return matriceRepere(vecteurZ, vecteurX, p);
}

Matrice Matrice::matriceRepere(Vecteur vecteurZ, Vecteur vecteurX, Point &pos)
{
	vecteurZ.normer();
	vecteurX.normer();
	Vecteur vecteurY = vecteurZ.pVectoriel(vecteurX);
	vecteurY.normer();
	vecteurX = vecteurY.pVectoriel(vecteurZ);
	vecteurX.normer();
	vecteurY = vecteurZ.pVectoriel(vecteurX);
	vecteurY.normer();

	return Matrice(vecteurX, vecteurY, vecteurZ, pos);
}